/**
 * Computer Graphics Lecture WS 2009/2010 Ulm University
 * Creator: Manuel Finckh
 * Email:   manuel.finckh@uni-ulm.de
 */

#include "render.h"

float ray_prob = 2.0;

float n1 = 1.0;
float n2 = 1.33;

float actual_n = n1;

Render::Render(Scene *scene) : scene(scene) {
  accel = new BVH(scene->triangles, scene->num_tris);

  cam = scene->cam;
  ResX = cam->ResX;
  ResY = cam->ResY;
  
  sResX = SNAPX;
  sResY = SNAPY;
  
  accum_buffer = new int[ResX * ResY];
  
  nr_samples = 0;

  image = new Vec3[ResX * ResY];
  for ( int i = 0; i < ResX * ResY; i++ ) image[i] = Vec3(0.0f, 0.0f, 0.0f);
  
  for ( int i = 0; i < ResX * ResY; i++ ) accum_buffer[i] = 0;
  
  snapshot = new Vec3[sResX * sResY];
  for ( int i = 0; i < sResX * sResY; i++ ) snapshot[i] = Vec3(0.0f, 0.0f, 0.0f);

  mtrand = new MTRand(1337);
  accum_index = 0;
}

Render::~Render() {
  delete accel;
  delete [] image;
  delete mtrand;
}

inline void Render::clearImage(){
	for ( int i = 0; i < ResX * ResY; i++ ){
		image[i] = Vec3(0.0f, 0.0f, 0.0f);
		accum_buffer[i] = 0;
	}
	nr_samples = 0;
}



float* Render::render(int shader, int snap) {
  float inv_accum, shrink;
  shrink_accum(inv_accum, shrink);
	
		actual = image;
		nResX = ResX;
		nResY = ResY;
	if(cam->moved){
		for ( int i = 0; i < nResX * nResY; i++ ) 
			actual[i] = Vec3(0.0f, 0.0f, 0.0f);
			clearImage();
	}
	depth = 0;

	  for ( int y = 0; y < nResY; y++ ) {
	    for ( int x = 0; x < nResX; x++ ) {
	    	if(shader == 4){
	    		int id;
	    		Vec3 light_point = scene->sampleRandomLight(id, mtrand->randExc(),mtrand->randExc(),mtrand->randExc());
	    		Vec3 dir2(0.0);
	    		//cout<<id<<endl;
			dir2 = sample_hemisphere(scene->triangles[id].getNormal());
			dir2.normalize();
			Ray ray(light_point,dir2,RAY_EPS, RAY_MAX);
		      	HitRec rec = accel->intersect(ray);
		      	// * scene->num_lights * scene->triangles[id].area()
		      	light_trace(ray,rec, scene->material[scene->mat_index[id]].color_e * scene->num_lights * scene->triangles[id].area());	
	    		nr_samples++;
	    	
	    	}else{	
			actual[x + y * nResX] *= shrink;
		     	float xc = x + mtrand->randExc(), yc = y + mtrand->randExc();
		     		
		     		Ray ray  = cam->getRay(xc,yc);
		      		HitRec rec = accel->intersect(ray);
				Vec3 hp = cam->eye + rec.dist * ray.dir;
				//float x1, y1;
				//cam->getPoint(hp-cam->eye,x1,y1);
				//cout<<"x: ("<<x<<" , "<<x1<<")   y:("<<y<<" , "<<y1<<")"<<endl;
		      		Vec3 color(0.0f, 0.0f, 0.0f);
		      		//cout<<"x: "<<x<<"   y: "<<y<<endl;
		      		if ( rec.id != -1 ) {
					switch ( shader ) {
				  	case 0: color = shade_debug(ray, rec); break;
				  	case 1: color = shade_simple(ray, rec); break;
				  	case 2: color = shade_path(ray, rec); break;
				  	case 3: color = shade_pathDL(ray, rec); break;
				  	case 5: color = shade_path_SSS(ray, rec); break;
				  	default: color = shade_simple(ray, rec); break;
				  	}
				}
			      	else if ( scene->environment != NULL ) {
					color = scene->getEnvironment(ray.dir);
				//ray.dir.printVec3();
			      	}
			      	actual[x + y * nResX] += color * inv_accum;
			      	actual_n = n1;
			}
			}
		}
		
		

  return (float*)actual;
}

Vec3 Render::shade_debug(Ray &ray, HitRec &rec) {
  //Vec3 normal = scene->triangles[rec.id].getNormal();
  Vec3 normal = scene->getShadingNormal(ray, rec.id);
  normal.abs();
  if ( scene->material[scene->mat_index[rec.id]].tex != NULL ) {
    // TODO
  }
  return normal;
}

Vec3 Render::shade_simple(Ray &ray, HitRec &rec) {
  //Vec3 normal = scene->triangles[rec.id].getNormal();
  Vec3 normal = scene->getShadingNormal(ray, rec.id);
  float cos = fabsf(normal * ray.dir);
  Vec3 color = cos * scene->material[scene->mat_index[rec.id]].color_d;
  if ( scene->material[scene->mat_index[rec.id]].tex != NULL ) {
    float alpha, beta, gamma;
    scene->triangles[rec.id].getBarycentric(ray, alpha, beta); gamma = 1.0f - alpha - beta;
    Vec2 coords = gamma * scene->uv[rec.id * 3] + alpha * scene->uv[rec.id * 3 + 1] + beta * scene->uv[rec.id * 3 + 2];
    
    color = Vec3::product(color, scene->material[scene->mat_index[rec.id]].tex->getColor(coords));
  }
  return color;
}

Vec3 Render::shade_mirror(Ray &ray, HitRec &rec){
	Vec3 normal = scene->getShadingNormal(ray, rec.id);
	Vec3 mirror = scene->triangles[rec.id].reflect(ray, normal);
	Vec3 hitpoint = ray.origin + rec.dist * ray.dir;
	Ray reflect_ray(hitpoint, mirror, RAY_EPS, RAY_MAX);
	HitRec r2 = accel->intersect(reflect_ray);

	      	Vec3 color(0.0f, 0.0f, 0.0f);
	      	if ( r2.id != -1 ) {
			return shade_mirror(reflect_ray, r2);
	      	}
	      	else if ( scene->environment != NULL ) {
			color = scene->getEnvironment(reflect_ray.dir);
		}	
		return color;
	
}

Vec3 Render::shade_path(Ray &ray, HitRec &rec){

	Vec3 color_e(0.0f);
	rec.setDepth(rec.depth+1);
	if(rec.id != -1 && rec.depth <= 5){
		
	
	
	
	/*if(scene->material[scene->mat_index[rec.id]].name.compare("Glass") == 0){
		//cout<<"test"<<endl;
		scene->material[scene->mat_index[rec.id]].glass(dir2, ray.dir, normal,mtrand->randExc());
	}
	else*/
	//if(mtrand->randExc() < ray_prob){
		//Vec3 normal = scene->triangles[rec.id].getNormal();
		Vec3 normal = scene->getShadingNormal(ray, rec.id);
		if(scene->material[scene->mat_index[rec.id]].color_e.length_sqr() > 0.0001){
			if(ray.dir * normal < 0.0)
				color_e = scene->material[scene->mat_index[rec.id]].color_e;
			else
				normal = -normal;
				
		}	
		//Vec3 normal = scene->getShadingNormal(ray, rec.id);
		
		//normal.normalize();
		Vec3 color = scene->material[scene->mat_index[rec.id]].color_d;
		//Vec3 mirror = scene->triangles[rec.id].reflect(ray, normal);
		Vec3 hitpoint = ray.origin + rec.dist * ray.dir;
	
	
	
		Vec3 dir2(0.0);
		dir2 = sample_hemisphere(normal);

		dir2.normalize();
		Ray reflect_ray(hitpoint, dir2, RAY_EPS, RAY_MAX);
		HitRec r2 = accel->intersect(reflect_ray);
		r2.setDepth(rec.depth);
		// texture ?
		/*if ( scene->material[scene->mat_index[rec.id]].tex != NULL ) {
	    		float alpha, beta, gamma;
	    		scene->triangles[rec.id].getBarycentric(ray, alpha, beta); gamma = 1.0f - alpha - beta;
	    		Vec2 coords = gamma * scene->uv[rec.id * 3] + alpha * scene->uv[rec.id * 3 + 1] + beta * scene->uv[rec.id * 3 + 2];
	    
	    		color = Vec3::product(color, scene->material[scene->mat_index[rec.id]].tex->getColor(coords));
	  	}*/
		      		
		      	// fabs(reflect_ray.dir*normal)*		
		return (color_e + Vec3::product(color, shade_path(reflect_ray, r2)));// / ray_prob;
	      	
	//}
	/*else if ( scene->environment != NULL ) {
		//color = scene->getEnvironment(ray.dir);
		color = Vec3(0.0);
	}	*/
	}
	return Vec3(0.0);

	
}

inline Vec3 Render::direct_illumination(Ray &ray, HitRec &rec){
	Vec3 hitpoint = ray.origin + rec.dist * ray.dir;
	int id;
	Vec3 light_dir = (scene->sampleRandomLight(id, mtrand->randExc(),mtrand->randExc(),mtrand->randExc()) - hitpoint);
	float dist = light_dir.length_sqr();
	light_dir.normalize();
	Vec3 light_normal = scene->triangles[id].getNormal();
	if(light_dir * light_normal < 0.0){
			Ray light_ray(hitpoint, light_dir, RAY_EPS, RAY_MAX);
			HitRec light_rec = accel->intersect(light_ray);
			
			// not in shadow
			if(light_rec.id == id){
				Vec3 radiance(0.0);
				//Vec3 normal = scene->triangles[rec.id].getNormal();
				Vec3 normal = scene->getShadingNormal(ray, rec.id);
				float geo_term = fabs(normal * light_dir) * fabs(light_normal * light_dir);
				geo_term /= dist;
				//cout<<light_dir.length_sqr()<<endl;
				if(scene->material[scene->mat_index[rec.id]].name.compare("SSS") == 0){
					radiance = scene->material[scene->mat_index[id]].color_e *  geo_term * 
					WardBRDF(-ray.dir,light_dir,normal)* scene->num_lights * scene->triangles[id].area();
				}else{
					radiance = scene->material[scene->mat_index[id]].color_e *  geo_term * 
					BRDF(-ray.dir,light_dir,normal)* scene->num_lights * scene->triangles[id].area();
				}
				//** scene->num_lights * scene->triangles[id].area() ;
				//cout<<geo_term<<endl;
				//radiance.printVec3();
				//scene->material[scene->mat_index[id]].color_e.printVec3();
				return radiance;
			}
		}

	return Vec3(0.0);


} 
// point + triangle id
inline Vec3 Render::direct_illumination(Vec3 & p, int &p_id){
	Vec3 hitpoint = p;
	int id;
	Vec3 light_dir = (scene->sampleRandomLight(id, mtrand->randExc(),mtrand->randExc(),mtrand->randExc()) - hitpoint);
	float dist = light_dir.length_sqr();
	light_dir.normalize();
	Vec3 light_normal = scene->triangles[id].getNormal();
	if(light_dir * light_normal < 0.0){
			Ray light_ray(hitpoint, light_dir, RAY_EPS, RAY_MAX);
			HitRec light_rec = accel->intersect(light_ray);
			
			// not in shadow
			if(light_rec.id == id){
				Vec3 radiance(0.0);
				//Vec3 normal = scene->triangles[p_id].getNormal();
				Vec3 normal = scene->getShadingNormal(light_ray, p_id);
				float geo_term = fabs(normal * light_dir) * fabs(light_normal * light_dir);
				geo_term /= dist;
				//cout<<light_dir.length_sqr()<<endl;
				radiance = scene->material[scene->mat_index[id]].color_e *  geo_term *(1/M_PI)* scene->num_lights * scene->triangles[id].area();
				return radiance;
			}
		}

	return Vec3(0.0);


}

inline Vec3 Render::shade_pathDL(Ray &ray, HitRec &rec){
	//cout<<rec.depth<<endl;
	Vec3 color(0.0f);
	rec.setDepth(rec.depth+1);
	if(rec.id != -1 && rec.depth <= 10){
		if(scene->material[scene->mat_index[rec.id]].color_e.length_sqr() > 0.0001){
			if(rec.depth <= 1 && (scene->triangles[rec.id].getNormal() * ray.dir) < 0.0)
				return scene->material[scene->mat_index[rec.id]].color_e;
			else
				return Vec3(0.0);
		}		
		//Vec3 normal = scene->triangles[rec.id].getNormal();
		Vec3 normal = scene->getShadingNormal(ray, rec.id);
		//color = scene->material[scene->mat_index[rec.id]].color_d;
		//color = Vec3(1.0);		

		Vec3 hitpoint = ray.origin + rec.dist * ray.dir;

		Vec3 direct = direct_illumination(ray, rec);
		
		//return Vec3::product(color,fabs(ray.dir*normal));	
		if(mtrand->randExc() < ray_prob){
			Vec3 dir2(0.0);
			color = scene->material[scene->mat_index[rec.id]].color_d;  
			
			if(scene->material[scene->mat_index[rec.id]].name.compare("Glass") == 0){
			//cout<<"glas"<<endl;
				float R, T;
				//cout<<"ray * normal: "<<ray.dir * normal<<endl;
				if(ray.dir * normal > 0.0){
					fresnel(R,T,ray.dir,-normal,n1,n2);
					//cout<<"von innen"<<endl;
					
				}
				else{
					fresnel(R,T,ray.dir,normal,n2,n1);
					//cout<<"von außen"<<endl;
				}
				//cout<<"R: "<<R<<"   T: "<<T<<endl;
				float r = mtrand->rand();
				//cout<<"r: "<<r<<endl;
				if( r < R)
					dir2 = scene->triangles[rec.id].reflect(ray);
						
				else{
					//cout<<"ref: ";
					if(ray.dir * normal > 0.0)
						refract(dir2, ray.dir, -normal, n1/n2);
					else
						refract(dir2, ray.dir, normal, n2/n1);
						
					//cout<<dir2*normal<<endl;;
				}
				direct = Vec3(0.0);
				color = Vec3(1.0);
				
			}
			else{
				//cout<<"kein glas"<<endl;
				dir2 = sample_hemisphere(normal);
			}
			//dir2 = sample_hemisphere(normal);
			dir2.normalize();
			Ray reflect_ray(hitpoint, dir2, RAY_EPS*10, RAY_MAX);
			HitRec r2 = accel->intersect(reflect_ray);
			r2.setDepth(rec.depth);
			      		
			
			//color.printVec3();
			//direct.printVec3(); 			
			return Vec3::product(color, (direct + shade_pathDL(reflect_ray, r2)) / ray_prob);
		}
		
		      	
	}
	else if ( scene->environment != NULL )
		//color = scene->getEnvironment(ray.dir);
		color = Vec3(0.0);
		
	return color;
}

inline float Render::rad2deg(float rad){

	return rad * 180.0 / M_PI;
}

inline void Render::light_trace(Ray &ray, HitRec &rec, Vec3 radiance){

	Vec3 color_e(0.0f);
	if(rec.id != -1){
		Vec3 normal = scene->triangles[rec.id].getNormal();
		if(scene->material[scene->mat_index[rec.id]].color_e.length_sqr() > 0.0001){
			if(ray.dir * normal < 0.0)
				color_e = scene->material[scene->mat_index[rec.id]].color_e;
			else
				normal = -normal;	
		}
		Vec3 hitpoint = ray.origin + rec.dist * ray.dir;
		Vec3 cam_dir = hitpoint - cam->eye;
		Ray shadow_ray(cam->eye, cam_dir, RAY_EPS, RAY_MAX);
	
		HitRec shadow_rec = accel->intersect(shadow_ray);
	
		// not in shadow
		
		if(shadow_rec.id == rec.id){
			float _x, _y, _weight;
			cam->getPoint(cam_dir, _x, _y, _weight);
	
			//cout<<_x<<"    "<<_y<<endl;
			//cout<<_weight<<endl;
		
			if(_x >= 0.0 && _x < cam->ResX && _y >= 0.0 && _y < cam->ResY){
				if(cam->w * cam_dir > 0.0){
					Vec3 color = scene->material[scene->mat_index[rec.id]].color_d;
					float dist = cam_dir.length_sqr();
					cam_dir.normalize();
					float geo_term = fabs(normal * cam_dir) * fabs(cam_dir * cam->w);
						geo_term /= dist;
				
					accumulate((int) _x, (int) _y,  Vec3::product((color), radiance * geo_term * BRDF(-ray.dir,cam_dir,normal) * _weight ));
				}	
					
			}
	
	
		}
		//if(mtrand->randExc() < ray_prob){
	
			Vec3 dir2(0.0);
			Vec3 color = scene->material[scene->mat_index[rec.id]].color_d;
			dir2 = sample_hemisphere(normal);
			dir2.normalize();
			Ray reflect_ray(hitpoint, dir2, RAY_EPS, RAY_MAX);
			HitRec r2 = accel->intersect(reflect_ray);
			float cos_l = fabs(dir2 * normal);
			light_trace(reflect_ray, r2, (Vec3::product(color, radiance) ));		
		//}
	}

}

inline Vec3 Render::shade_path_iterative(Ray &ray, HitRec &rec) {
Vec3 color = Vec3(0.0f);
  Vec3 path_weight(1.0f);
  int depth = 0;
  
  while ( true ) {
    Vec3 point = ray.origin + rec.dist * ray.dir;
    Vec3 s_normal = scene->getShadingNormal(ray, rec.id);
    Vec3 g_normal = scene->triangles[rec.id].getNormal();
    Material *mat = &scene->material[scene->mat_index[rec.id]];

    float cos_g = g_normal * ray.dir;
    if ( cos_g > 0.0f ) g_normal = -g_normal;
    if ( s_normal * g_normal < 0.0f ) s_normal = -s_normal;
    float cos_s = -(s_normal * ray.dir);
    if ( cos_s < 0.0f ) break;
    if ( cos_g < 0.0f) {
      color += Vec3::product(mat->color_e, path_weight);
    }
    if ( ++depth > 5 ) break;

    path_weight = Vec3::product(path_weight, mat->color_d);
   
    Vec3 dir_new;
    Material::diffuse(dir_new, s_normal, mtrand->rand(), mtrand->rand());
    ray.origin = point;
    ray.dir = dir_new;
    ray.tmin = RAY_EPS;
    ray.tmax = FLT_MAX;

    rec = accel->intersect(ray);
    if (rec.id == -1 ) break;
  }
  if ( scene->environment != NULL ) {
    color += Vec3::product(path_weight, scene->getEnvironment(ray.dir));
  }
  return color;
}

inline float Render::rand_exponential(const float &lambda){
	return (-1.0/lambda) * log(mtrand->rand());
}

inline float Render::schlick(const Vec3 &in, const Vec3 &n){


	float theta = fabs(in * n);
	float r0 = ((n1 - n2) / (n1 + n2)) * ((n1 - n2) / (n1 + n2));
	
	return (r0 + (1.0-r0) * pow((1.0-cos(theta)),5.0));
}

float phase = 1.0 / (4 * M_PI);

float albedo = 0.5;

float sigma_s = 2.0;
float sigma_a = 0.00035;
float sigma_t = sigma_s + sigma_a;
float sigma_tr = sqrtf(3.0*sigma_a * sigma_t);

inline float Render::dipole(Vec3 x0, Vec3 normal , Vec3 &x1, int &id){

	float Fdr = 1.440 / (n2*n2) + 0.710 / n2 + 0.668 + 0.0636*n2; 
	
	float albedo_r = sigma_s / sigma_t;
	
	float zr = 1.0 / sigma_t;
	float A = (1.0 + Fdr) / (1.0 - Fdr);
	float D = 1.0 / (3.0 * sigma_t);
	float zv = zr + 4.0 * A * D;
	
	Vec3 u,v,w;
	Vec3::onb(u,v,w,normal);
	int z = 1;
	HitRec v1;
	Ray upper;
	float d;
	do{
		d = rand_exponential(sigma_tr);
		//cout<<d<<endl;
		Vec3 l = mtrand->rand() * u + mtrand->rand() * v;
		l.normalize();
		x1 = x0 + d * l;
	
		upper = Ray(x1 + z * normal, -normal, RAY_EPS, RAY_MAX);
		//Ray lower(x0 -  normal, x1 - (x0 -  normal), RAY_EPS, RAY_MAX);
		v1 = accel->intersect(upper);
	}while(v1.id == -1);
	
	id = v1.id;
	x1 = upper.origin + v1.dist * upper.dir;
	d = (x1-x0).length();	
	//cout<<d<<endl;
	
	float dr = sqrtf(zr*zr + d*d); 
	float dv = sqrtf(zv*zv + d*d);
	float term1 = (sigma_tr * dr + 1) * exp(-sigma_tr*dr)/(sigma_t * pow(dr,3.0));
	float term2 = zv*(sigma_tr * dv + 1) * exp(-sigma_tr*dv)/(sigma_t * pow(dv,3.0));
	
	return (albedo_r / (4.0 * M_PI) * (term1 + term2));
	
	

}


inline Vec3 Render::shade_path_SSS(Ray &ray, HitRec &rec){
	//cout<<rec.depth<<endl;
	Vec3 color(0.0f);
	rec.setDepth(rec.depth+1);
	if(rec.id != -1 && rec.depth <= 5){
		if(scene->material[scene->mat_index[rec.id]].color_e.length_sqr() > 0.0001){
			if(rec.depth <= 1 && (scene->triangles[rec.id].getNormal() * ray.dir) < 0.0)
				return scene->material[scene->mat_index[rec.id]].color_e;
			else
				return Vec3(0.0);
		}		
		//Vec3 normal = scene->triangles[rec.id].getNormal();
		Vec3 normal = scene->getShadingNormal(ray, rec.id);
		

		Vec3 hitpoint = ray.origin + rec.dist * ray.dir;
		Vec3 SSS(0.0);
	
		
		if(scene->material[scene->mat_index[rec.id]].name.compare("SSS") == 0){
			int z = 5;
			
			float s0 = rand_exponential(sigma_t);
			Vec3 refr(0.0);
			
			HitRec inner;
			Ray bounce_ray;
			refract(bounce_ray.dir,ray.dir,normal, n1/n2);
			inner.dist = s0;
			bounce_ray.origin = hitpoint;
			while(z--){
				Vec3 P = bounce_ray.origin + inner.dist * bounce_ray.dir;
				int id;
				Vec3 light_P = scene->sampleRandomLight(id, mtrand->randExc(),mtrand->randExc(),mtrand->randExc());
				
				Ray light_ray(light_P, P - light_P, RAY_EPS, RAY_MAX);
				HitRec light_rec = accel->intersect(light_ray);
			
			
				if(light_rec.id >= scene->material[scene->mat_index[rec.id]].start && light_rec.id <= scene->material[scene->mat_index[rec.id]].end){
					//normal = scene->triangles[light_rec.id].getNormal();
					normal = scene->getShadingNormal(ray, light_rec.id);
					Vec3 in_P = light_ray.origin + light_rec.dist * light_ray.dir;
					float s1 = (in_P - P).length();
					// nachher
					float c_i = fabs(light_ray.dir * normal);
					float s1_prime = s1 * (c_i / (sqrtf(1.0 - (1.0/n2)*(1.0/n2) * (1.0 - (c_i * c_i)))));
				
					//Vec3 direct(0.0);
				
					float t1, t2, r;
					fresnel(r,t1,ray.dir,normal,n2,n1);
					fresnel(r,t2,light_ray.dir,normal,n2,n1);
				
					//float geo_term = 1.0; //fabs(normal * ray.dir) / fabs(normal * light_ray.dir);
					Vec3 light_normal = scene->triangles[id].getNormal();
					float geo_term = fabs(light_normal * light_ray.dir);
					geo_term /= (light_rec.dist * light_rec.dist);
					float sigma_tc = sigma_t;// + (fabs(normal * ray.dir) / fabs(normal * light_ray.dir)) * sigma_t;
					//cout<<light_dir.length_sqr()<<endl;
				
					SSS = scene->material[scene->mat_index[id]].color_d * scene->num_lights * geo_term * scene->triangles[id].area();
					SSS *= (sigma_s * t1 * t2 * phase / sigma_tc) * exp(-s1_prime*sigma_t) * exp(-s0*sigma_t) / M_PI; 
					color += SSS;
					s0 = rand_exponential(sigma_t);
				
				
					bounce_ray = Ray(P,sample_sphere(), 0.0, RAY_MAX);
					inner = accel->intersect(bounce_ray);
					if(inner.id != -1){
						if(inner.dist < s0){
							z = 0;
							// später vllt: refract
							hitpoint = bounce_ray.origin + inner.dist * bounce_ray.dir;
							normal = scene->triangles[inner.id].getNormal();
							rec = inner;
							ray = bounce_ray;
						}else{
							inner.dist = s0;
						}
				
					}	
				}
		     	}		
		} 
		
		Vec3 dir2(0.0);
			dir2 = sample_hemisphere(normal);
			dir2.normalize();
		// get direct illumination
		Vec3 direct = direct_illumination(ray, rec);
		/*
		// dipol zeug
		if(scene->material[scene->mat_index[rec.id]].name.compare("SSS") == 0){
		float t1, t2,r;
		//fresnel(r,t1,ray.dir,normal,n1,n2);
		//fresnel(r,t2,dir2,normal,n1,n2);
		t1 = t2 = 1.0;
		int x_id;
		Vec3 x1;
		float Rd = dipole(hitpoint,normal, x1, x_id);
		Vec3 x_direct = direct_illumination(hitpoint,x_id);
		direct += Rd * x_direct * t1 * t2 * (1.0/M_PI) * (1.0 / sigma_tr);
		}*/

	
		//if(mtrand->randExc() < ray_prob){
			
			Ray reflect_ray(hitpoint, dir2, RAY_EPS, RAY_MAX);
			HitRec r2 = accel->intersect(reflect_ray);
			r2.setDepth(rec.depth); 	
			color = scene->material[scene->mat_index[rec.id]].color_d;
			return Vec3::product(color, (direct  + shade_path_SSS(reflect_ray, r2)));
		//}
	
		      	
	}
	else if ( scene->environment != NULL )
		color = Vec3(0.0);
		
	return color;
}


// n1 = actual medium, n2 = adjacent one
inline void Render::fresnel(float &R, float &T, const Vec3 &in, const Vec3 &n, float n1, float n2){

	float n0 = n1/n2;
	float c_i = fabs(in * n);
	float s_t = n0 * n0 * (1.0 - c_i*c_i);
	float c_t = sqrtf(1.0 - s_t);
	//cout<<"critical: "<<rad2deg(asin(n0))<<"  winkel: "<<rad2deg(acos(c_i))<<endl;
	if(s_t > 1.0){
		R = 1.0;
		T = 0.0;
		return;
	}	
	
	float R_r = (n1 * c_i - n2 * c_t) / (n1 * c_i + n2 * c_t);
	float R_t = (n2 * c_i - n1 * c_t) / (n2 * c_i + n1 * c_t);
	
	R = 0.5 * (R_r * R_r + R_t * R_t);
	T = 1.0 - R; 
	


}


inline void Render::refract(Vec3 &out, const Vec3 &in, const Vec3 &n, float n0){
	

		float c1 = n * (-in);
		float c2 = sqrt(1 - (n0 * n0) * (1-(c1*c1)));
		//if(c1 > 0)
			out = n0 * in + (n0 * c1 - c2) * n;
		//else
		//	out = n0 * in + (n0 * c1 + c2) * n;
			
}

inline void Render::shrink_accum(float &inv_accum, float &shrink) {
  if (!cam->moved) {
    accum_index++;
    inv_accum = 1.0f / (float)accum_index;
    shrink = (accum_index-1) * inv_accum;
  } else {
    accum_index = 1;
    inv_accum = 1.0f;
    shrink = 0.0f;
    for(int i = 0; i<ResX*ResY; i++)
    	accum_buffer[i] = 1;
    clearImage();
    
  }
}

inline void Render::accumulate(const int &x, const int &y, const Vec3 &color){

	/*
	Vec3 t = actual[x + y * ResX];
	int N = accum_buffer[x + y * ResX];
	actual[x + y * ResX] = (N * t + color) / (N+1);
	accum_buffer[x + y * ResX]++;
	*/
	
	Vec3 t = actual[x + y * ResX];
	int N = nr_samples;
	actual[x + y * ResX] = (N * t + color) / (N+1);
	

}


inline Vec3 Render::sample_hemisphere(Vec3 normal){
	/*float theta = mtrand->randExc(2.0 * M_PI);
	float phi = mtrand->randExc(M_PI / 2.0);
	float r = cos(phi);
	Vec3 tmp(r*cos(theta),r*sin(theta),sin(phi));
	tmp.normalize();*/
	float u = mtrand->randExc();
	float v = mtrand->randExc();
	Vec3 u1, v1, w1, out;
  Vec3::onb(u1, v1, w1, normal);

  const float phi = 2.0f * M_PI * u;
  const float sqrtu = sqrtf(v);
  const float x = sqrtu * cosf(phi);
  const float y = sqrtu * sinf(phi);
  const float z = sqrtf(1.0f - v);

  out = x * u1 + y * v1 + z * w1;
	return out;
}

inline Vec3 Render::sample_sphere(){

	float u = mtrand->randExc();
	float v = mtrand->randExc();
	Vec3 out(0.0);


	const float phi = 2.0f * M_PI * u;
	const float theta = M_PI * v;
	out.x = sin(theta) * cos(phi);
	out.y = sin(theta) * sin(phi);
	out.z = cos(theta);
	return out;
}


float alpha = 0.071, rho_s = 0.07, rho_d = 0.60;
float alpha_s = alpha * alpha;
// Ward BRDF
inline float Render::WardBRDF(Vec3 in, Vec3 out, Vec3 normal){
	
	in.normalize();
	out.normalize();
	Vec3 half = in + out;
	half.normalize();
	Vec3 x1, y1, z1;
  	Vec3::onb(x1, y1, z1, normal);
	float gamma = tan(acos(half * normal));
	float s = rho_s / (4.0 * alpha_s * M_PI * sqrtf(fabs( (in * normal) * (out * normal))));
//std::cout<<"BRDF(s): "<<4.0 * alpha * alpha *  (in * normal) * (out * normal)<<std::endl;
	//float exponent = ((half * x1) * (half * x1) / alpha_s + (half * y1) * (half * y1)/ alpha_s)  / ((half * normal)*(half * normal));
	float exponent = gamma * gamma / alpha_s;
	//std::cout<<"BRDF: "<<((rho_d + s * exp(-exponent)))<<std::endl;
	return (rho_d/M_PI + s * exp(-exponent));
	//return 1.0/M_PI;
	//return 1.0;

}

inline float Render::BRDF(Vec3 in, Vec3 out, Vec3 normal){
	
	return 1.0/M_PI;

}





